function [dof_map, V, T, u1, u2] = Lshape_Crack(n_round, Young, nu, mode, FE_Order, method, FE_Type)
% function [dof_map, V, T, u1, u2] = Lshape_Crack(n_round, Young, nu, mode, FE_Order, method, FE_Type)
%
%  Solve planar elasitic problem on domain L: [-1,1]x[-1,1]/[0,1]x[0,1] 
%  with analytic solution given by fun_u
%  the right hans side is then given by fun_f
%
%  see the example in paper of M. Ainsworth and B. Senior:
%  Comput. Methods Appl. Mech. Engrg. 150(1997), 65-87.
%
%  This is also for validation!!
%
%  Test pass on Aug. 11th, 2012
%
%  Author: Dr. Xian-Liang Hu
%

if nargin < 2
    Young = 1;   % Young's modulus (N/mm^2)
    nu = 0.3;               % Poisson's coefficient , whose limit is 0.5
end

if nargin < 4
    FE_Order = 3;   % default polynomial order.
end

if nargin < 5
    method = 2;  % using FEM style by default
end

if nargin < 6
    FE_Type = 'BB';  % default FE basis using Bernstein Polynomial
end


% %% for output dofs v.s. errors
% dof_count = zeros(10,1);
% err_L2 = zeros(10,1);
% err_inf = zeros(10,1);

Quad_Order = 13;

% 1: prepare the mesh
[p, ~, t] = initmesh('lshapeg','Hmax',0.35,'Hgrad',1.99);
V = [p(1,:)' -p(2,:)']; T = t(1:3,:)';
[T, E, ET, TE] = build_fem_mesh(V, T);
[V, T, TE, E, ET] = refine_L_shape(V, T, TE, E, ET, 8);

V = V*sqrt(2);  % do rescale
x = V(:,1); y = V(:,2); theta = 45/180*pi;  
V(:,1) =  cos(theta)*x + sin(theta)*y;  % do rotate
V(:,2) = -sin(theta)*x + cos(theta)*y;

for round = 1:n_round  % do global refinement
    [V,T, TE, E, ET] = refine_mesh_uniform(V, T, TE, E, ET);
end

fprintf('There are %d triangles.\n', size(T,1));


%%%%%%%%
% 2: simply posing Dirichlet condition for all the boundaries
% bdr_Neumann, @fun_f are actually dummy
% 
bdr_Dirichlet = find(ET(:,2)==0);
bdr_Neumann = [];

%%%%%%%%%%%%%
% 3: calling FEM routine for two dimensional elasitic problems
[dof_map, u1, u2, t] = elastic_fem(V, T, TE, ET, FE_Type, FE_Order, Quad_Order, method, ...
                                     Young, nu, @fun_f, bdr_Dirichlet, @fun_u, bdr_Neumann, @fun_surface_f, mode);

                                 
%%% Trace the time Consuming
fprintf('Totally cost %5.4f seconds.\n', t(end) - t(1));
                                 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  The right hand side functions are ZERO in these cases.
%
% 
function [f_1, f_2] = fun_f(xx, yy, E, nu, varargin)
    f_1 = E*zeros(size(xx));
    f_2 = nu*zeros(size(yy));
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  The surface force function applied on Neumann boundary
% 
% !!! Important !!! Note that this is a special example:
% Since the exact solution is known, we can simply impose the Dirichlet
% boundary condition on all the boundary. So that the surface force is set
% to be ZERO here!
%
function [fb_1, fb_2] = fun_surface_f(xx, yy, E, nu, varargin)
    fb_1 = E*zeros(size(xx));
    fb_2 = nu*zeros(size(yy));
end